Heat conduction in graphene flakes with inhomogeneous mass interface 
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Using nonequilibrium molecular dynamics simulations, we study the heat conduction in graphene flakes 
composed by two regions. One region is mass-loaded and the other one is intact. It is found that the mass 
interface between the two regions greatly decreases the thermal conductivity, but it would not bring thermal 
rectification effect. The dependence of thermal conductivity upon the heat flux and the mass difference ratio 
are studied to confirm the generality of the result. The interfacial scattering of solitons is studied to explain the 
absence of rectification effect. 
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Thermal rectification is a phenomenon that the heat flux 
runs preferentially along one direction and inferiorly along the 
opposite direction[L- 14|. It has attracted a great deal of atten- 
tion in the last decade since it reveals the possibility to control 
the heat transportation process. With an improved understand- 
ing of thermal rectification, various devices like thermal tran- 
sistors, thermal logic circuits and thermal diodes could be fab- 
ricated. Two methods are commonly used to design thermal 
rectifiers. The first method is to couple two or more anhar- 
monic chains with different nonlinear potentials together|3]- 
Q. The explanation for the observed rectification effect is that 
the phonon bands of different regions of the chain will change 
from overlap to separation when the heat flux is reversed. The 
asymmetry of interaction potential controls the phonon band 
shift and it plays the central role here. The second method 
is to implement asymmetric geometric shape in quasi- ID and 
2D systems. For example, it is applied in deformed carbon 
nanotubes|6|, carbon nanohoms|7 1, triangle shaped, trapezoid 
shaped and U-shaped graphene flakeslS nlOJ . Thermal con- 
ductivity is higher when the heat flux runs from the narrow to 
the wide region. The explanation for the observed rectification 
is that the asymmetric geometric shape introduces asymmetric 
boundary scattering of phonons. The asymmetry of geometric 
shape controls the phonon scattering and it plays the central 
role in this case. 

Recently a new procedure is considered by Chang et al. in 
carbon and boron nitride nanotubesfTTl. They introduce the 
asymmetry of mass distribution by covering external platinum 
compound on the nanotubes. Comparing with nanotubes, the 
thermal contribution of the platinum compound can be ne- 
glected. So the mass-loading procedure is idealized by Chang 
et al. as changing the atomic mass of the atoms in the heat 
conduction process. Higher thermal conductivity is observed 
when the heat flux runs from the heavy mass region to the 
light mass region. However the observed rectification effect 
cannot be explained by phonons. First, the externally loaded 
molecules do not contribute to any bond, thus the anharmonic 
coupling between the atoms in the nanotubes is not changed. 
The asymmetry of interaction potential does not contribute 
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here. Second, the associated geometric deformation of the 
nanotubes is not dominant, thus the boundary scattering of 
phonons is also not changed. The asymmetry of geometric 
shape also does not contribute in this case. Then Chang et 
al. surmise that the origin of the rectification effect is due to 
the asymmetric interfacial scattering between the two regions 
with inhomogeneous mass. However as pointed by Change et 
al, it is well known that the interfacial scattering of phonons 
is independent on the incident direction. Therefore they pos- 
tulate that solitons are responsible for the asymmetric inter- 
facial scattering process. Later similar thermal rectification 
effect has been observed in molecular dynamic simulation of 
mass-graded carbon nanotubes fT2'-T4l . In the simulations, 
the atomic mass of carbon atoms is gradually increased from 
12 to 300 along the axis of tubes. The setup is treated as a 
combination of multiple inhomogeneous mass interfaces. 

Graphene|15, 16|, a single layer of carbon atoms in a hon- 
eycomb lattice with sp^ bonds, reveals superior high thermal 
conductivity up to 2500-5000 W/mK at room temperature 1 17] 
[TSl . It has been considered as a promising candidate for var- 
ious thermal devices. Graphene is similar to carbon nanotube 
in both structure and thermal property, so it is interesting to 
know whether thermal rectification would occur if the mass 
interface is implemented. Besides, managing heat transporta- 
tion across the interface is very common in nanoelectronic de- 
sign, thus it is also interesting to understand the influence of 
mass interface upon the thermal conductivity of graphene. 

Here we study the heat conduction in graphene flakes 
with inhomogeneous mass interface between two regions by 
NEMD (nonequilibrium molecular dynamics) simulations. 
The graphene flakes are composed by two regions. One re- 
gion is intact and the other region is externally mass-loaded. 
For simplicity the mass-loaded region is realized by modulat- 
ing the atomic mass of carbon atoms in the molecular dynam- 
ics simulations. This setup is widely accepted as the ideal- 
ized method for the mass-loading procedure 1, 1 2-,- 1 4 J . We re- 
port that thermal rectification effect do not exist even a large 
mass difference ratio is implemented. In order to understand 
the microscopic mechanism leading to the absence of thermal 
rectification effect, we also study the interfacial scattering of 
solitons in graphene. 

We carry out the molecular dynamics simulations of the 
heat conduction in two rectangle graphene flakes with zig-zag 
edges along the x-axis and armchair edges along the y-axis. 
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Figure 1 : (a) Schematic of the graphene flakes with inhomogeneous 
mass interface. The first one is labeled as the niQ — niff graphene 
flake. The upper region of the mo — niff graphene flake is mass- 
loaded which is simplified by modulating the atomic mass of the 
carbon atoms as 12. Similarly the second one is labeled as the 

m^r — niQ graphene flake. The lower region of the mif — mo graphene 
flake is mass-loaded, (b) The typical temperature profiles of the two 
graphene flakes in the heat conduction process. Here mf^ /mo = 4 and 
J = 0.35 eV/ps are implemented. The 12 — 48 and 48—12 graphene 
flakes exhibit the same temperature drop near the mass interface and 
between the two ends along the whole simulation cell. 



The graphene flakes are shown in Fig. 1(a). They are both 
288 A long and 21 A wide. The heat sources and heat sinks 
are covered by red boxes with 100 carbon atoms. The outmost 
edges of the heat sources/sinks are frozen. It corresponds to 
fixed boundary conditions in the y-axis. Periodic boundary 
condition is used along the x-axis. The atomic mass of the 
intact carbon atoms is mo = 12 and they are drawn in cyan. 
The atomic mass of the mass-loaded carbon atoms is niN ^ 12 
and they are drawn in silver. The upper half region of the first 
graphene flake in Fig. 1(a) is mass-loaded, thus we label it as 
the as the mo — mm graphene flake. Similarly we label the sec- 
ond one in Fig. 1(a) as the myv — mo graphene flake. The heat 
flux runs along the m^r — mo graphene flake is equivalent to 
the reversed heat flux runs along the mo — m^ graphene flake. 
Thermal conduction process is investigated by imposing the 
same heat flux along the two graphene flakes. It is much more 
convenient later to compare the temperature profiles since the 
heat sources and heat sinks are in the same direction. We use 
the same reactive empirical bond-order (REBO) potential lfT9l 
as implemented in the LAMMPS 1 20 1 code to simulate the an- 
harmonic coupling between the carbon atoms. Equations of 
motions are integrated with velocity Verlet algorithm with the 
minimum timestep At — 0.25 fs. 

First the graphene flakes are equilibrated at a constant tem- 
perature T=300 K in the Nose-Hoover thermostat by 0.75 ns. 
After that the heat flux is imposed. It is realized by the energy 
and momentum conserving velocity rescaling algorithm de- 
veloped by Jude and Jullien ll2Ti . A constant rate of kinetic 
energy dE is added in the heat source and removed in the 
heat sink at each time step dt. The heat flux is calculated as 
J = dE /dt. It is widely used to study the interfacial heat con- 
duction in different materials [i22n24J . We divide the graphene 



flakes by several 4 A long slabs along the y-axis to obtain the 
temperature profiles. Each slab contains about 60-70 carbon 
atoms. The local temperature in each slab is calculated from 
the averaged kinetic energy of the carbon atoms. The tem- 
perature profiles are averaged over every 100 ps after the heat 
flux is imposed. The whole nonequilibrium simulation pro- 
cess covers 3 ns and the last 1 ns is utilized as the steady state 
since the temperature profiles do not change much. Thermal 
conductivity G of the whole graphene flake is calculated by 
the Fourier's Law: 
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Here J is the heat flux, A is the cross section of the heat trans- 
fer defined by the product of the width and the thickness of the 
graphene flakes. The thickness of the graphene flake is usu- 
ally considered as the bond length of the carbon atoms and it 
is taken as 1.4 A in our computationf^S^/SSl. AT (AL) is the 
temperature (distance) difference between the two ends on the 
graphene flake. Thermal conductivity G represents the ther- 
mal conducting capacity of the whole graphene flake. eV/ps 
is selected as the unit for the heat flux, A is selected as the 
unit for the width and thickness, K is selected as the unit for 
the temperature. W/mK is selected as the unit for the thermal 
conductivity, so the associated numerical results are converted 
to this unit. 

In Fig. 1(b) we show the typical temperature profiles of the 
two graphene flakes. Here m^/mo — 4 and J ^ 0.35 eV/ps are 
implemented. A sudden temperature drop is observed near the 
mass inteif ace. It indicates that the mass interface behaves like 
a strong thermal resistive boundary. In Fig. 1(b), it also shows 
that although the temperature profiles are different in the two 
graphene flakes, but the temperature differences between the 
two ends are the same. The temperature difference in the 
12 — 48 graphene flake is 94.0 K and the temperature differ- 
ence in the 48 — 12 graphene flake is 93.8 K. Thus the thermal 
conductivity G of the two graphene flakes are the same. It in- 
dicates that even for such a large temperature difference, there 
is still no observable thermal rectification effect. The 12 — 48 
and the 48 — 12 graphene flakes have the same thermal con- 
ducting capacity and the heat flux runs equivalently without 
preferred direction. 

In order to further confirm the result that there is no ther- 
mal rectification effect brought by the mass interface in the 
graphene flakes, the dependence of thermal conductivity upon 
the heat flux and the mass difference ratio are studied. First we 
keep the mass difference ratio m^ /mo — 4 invaiiant and vary 
the heat flux from 0.15 to 0.5 eV/ps. In Fig. 2(a) we show 
the dependence of thermal conductivity upon the heat flux. 
Thermal conductivity is almost the same in the two graphene 
flakes. It indicates that there is no thermal rectification effect 
by using different value of heat flux. Meanwhile it is found 
that thermal conductivity is decreasing with the heat flux. For 
/ = 0.15 eV/ps, the thermal conductivity is about 1 10 W/mK. 
For / = 0.5 eV/ps, it is reduced to 38 W/mK which is about 
35% of the original value. It indicates that interfacial scat- 
tering is enhanced by the amount of heat flux and the effect 
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Figure 2: (a) Thermal conductivity G vs heat flux J. Here the mass 
difference ratio inj\i /niQ = 4 is unchanged and the heat flux / varies 
from 0.15 to 0.5 eV/ps. (b) Thermal conductivity G vs the mass dif- 
ference ratio mif /mo. Here the heat flux i = 0.35 eV/ps is unchanged 
and the mass difference ratio varies from 1 to 5. For mt^/niQ = 1, the 
thermal conductivity G is obtained from the graphene flake without 
the mass interface. 



might be taken into consideration in real application. 

Second we keep the heat flux J = 0.35 eV/ps invariant and 
vary the mass difference ratio 111^/1110 from 1 to 5. Here 
inN/mo = 1 stands for the graphene flake without the mass 
interface. In Fig. 2(b) we show the dependence of thermal 
conductivity upon the mass difference ratio. Thermal con- 
ductivity is almost the same in the two graphene flakes. It 
indicates there is no thermal rectification effect by using dif- 
ferent mass difference ratio. Meanwhile it is found that ther- 
mal conductivity is greatly decreased by the mass interface. 
The thermal conductivity of the graphene flake without the 
mass interface is about 269 W/mK. When the mass interface 
is implemented (111^/1110 = 5), it is reduced to 45 W/mK. It 
is only about 16% of the original value. Here it provides a 
possible route to tune the thermal behavior of graphene. The 
mass interface can be implemented by two methods. One is 
to load external heavy and thermal insulating molecules upon 
the carbon atomslfTT Hl4l . The other one is to use different 
ratio of isotope substitutions ll27l - E9l which is demonstrated 
possible in experiment by chemical vapor deposition growth 
of graphene on metal li30i . 

The heat conduction results suggest that inhomogeneous 
mass interface cannot lead to thermal rectification effect in 
graphene. In the above simulations, the anharmonic coupling 
between the carbon atoms is kept constant and there is no 
geometric deformation by using periodic boundary condition 
along the width. Therefore just like carbon nanotubes, the 
asymmetry of interaction potential and the asymmetry of ge- 
ometric shape also do not contribute in the graphene flakes. 
Furthermore, the interfacial scattering of solitons in surmised 
to be responsible for the thermal rectification effect in car- 
bon nanotubes [TT'l. So in order to understand the absence of 
thermal rectification effect, it is necessary to study the inter- 
facial scattering of solitons in graphene. Recently, subsonic 
NLS (Nonlinear Schrodinger) equation described solitons are 
found in graphenel31 1. They preserve their identities in prop- 
agation and exhibit strong interactions and phase shifts in col- 
lision. Their energy reflection rates between the two regions 
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Figure 3: [(a)-(b)] The barrier line in the middle of (a) and (b) sep- 
arates the graphene flake into two regions. The atomic mass of the 
mo region is 12. The atomic mass of the ra^v region is 48. The curve 
ti represents the longitudinal velocity of the carbon atoms before the 
scattering. The curve t2 represents the longitudinal velocity of the 
carbon atoms after the scattering. The time interval between ti and 
t2 is 3.2 ps. The 12 — 48 scattering is shown in (a). The 48 — 12 scat- 
tering is shown in (b). [(c)-(d)] Rescaling the width of the reflected 
soliton S2 and the transmitted soliton S3 according to the width of the 
incident soliton si . The rescaling in the 12 — 48 scattering is shown 
in (c). The rescaling in the 48 — 12 scattering is shown in (d). 



with different mass are needed if one tends to know whether 
the solitons would bring rectification effect in the interfacial 
scattering process ll32l434l . 

Here we first generate a longitudinal soliton in a 485 nm 
long graphene flake with atomic mass mo = 12. The chirality 
of the graphene flake is the same as the two graphene flakes 
in Fig. 1(a). In the curve ti of Fig. 3(a), we show the longitu- 
dinal velocity (vy) distribution of the carbon atoms before the 
interfacial scattering. The soliton is set as the incident soliton 
and its amplitude is positive. After that we change the atomic 
mass of the carbon atoms ahead of the solitons to be = 48 
to set up a mass interface. The barrier line in the middle of Fig. 
3(a) separates the graphene flake into two regions. The atomic 
mass of the left region is intact (mo — 12) and the atomic mass 
of the right region is mass-loaded (;?% = 48). The interfacial 
scattering from the mo to the m^ region happens when the in- 
cident soliton hits the mass interface. In the curve t2 of Fig. 
3(a), we show the longitudinal velocity distribution of the car- 
bon atoms after the interfacial scattering. The incident soliton 
is scattered to a transmitted soliton and a reflected soliton. The 
amplitude of the transmitted soliton is positive and the ampli- 
tude of the reflected soliton is negative in the 12 — 48 scat- 
tering. Similarly we obtain the interfacial scattering from the 
mfi to the mo region. In Fig. 3(b) we show that the amplitude 
of the transmitted soliton is positive and the amplitude of the 
reflected soliton is also positive in the 48 — 12 scattering. The 
scattering results indicate that the amplitude of the reflected 
soliton is dependent upon the incident direction. It is negative 
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in the mo — scattering and positive in the mpi — mo scat- 
tering. Such behavior of the reflected soUton has also been 
obtained in the interfacial scattering of similar NLS solitons 
in ID nonlinear chain[32|. It represents a different kinetic be- 
havior of the NLS solitons from the KdV solitons which have 
no reflected solitons in the — mo scattering. 

There is a rescaling relation between the width of the scat- 
tered solitons. The propagating velocity of a soliton in the 
mo region is 20 km/s and in the m^ section llSTl is v{mN) — 
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tering time At is very short thus we can neglect the amplitude 
dependent parameters in the formula to estimate Af . In the 
mo — m^ scattering Af can be estimated by the width of the 
incident soliton w{,„ „, as Af 
the reflected soliton 



"10 - 

transmitted sohton wJ^Q-mN 



m,-mJ^Q- So the width of 

in the mo region and the width of 



in the m^ region are: 



w, 



v{mpi)/St - 



mo J 
w 



(2) 



(3) 



Similar relations can be obtained in m^ — mo scattering. The 
width of the reflected soliton w^^_,„g in the m^ region and the 
width of transmitted soliton in the mo region are: 
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(4) 
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In Fig. 3(c) and Fig. 3(d) we show the width rescaling rela- 
tions between the scattered solitons. The width of the reflected 
soliton S2 and the transmitted soliton S3 are rescaled well ac- 
cording to the width of the incident soliton s 1 . 

In order to understand the role of solitons play in the heat 
conduction process, the energy reflection rate is needed. We 
measure the energy E and the momentum P of a soliton in 
graphene as the aggregated momentum and energy of all the 
carbon atoms along the width of the soliton ll35ll36]l : 
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Here only the carbon atoms with || Vi || > 0.01 m/s are counted. 
The widths of the scattered solitons also could be estimated by 
the rescaling relation in Eq. (2)-(5). When the energy and the 
momentum of the scattered solitons are obtained, the energy 
reflection rate and the momentum reflection rate are 
defined as: 
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Figure 4: Here mo = 12 and m^r varies from 12 to 60. For mj\i/ino = 
1 , it stands for the case that there is no mass interface in the graphene 
flake. Thus the energy and the momentum reflection rates are both 
in this case, (a) Energy reflection rate vs mass difference ratio 
OTiv/mo. (b) Momentum reflection rate vs mass difference ratio 
m^r /mo . 



Here E' and P' are the energy and the momentum of the inci- 
dent soliton. E'^ and are the energy and the momentum of 
the reflected soUton. 

In Fig. 4 we show the dependence of the energy and the mo- 
mentum reflection rate upon the mass difference ratio. In Fig. 
4(a) it shows that the energy reflection rates in the mo — m^ 
and the mf/ — mo scattering are the same. In Fig. 4(b) it 
shows that the momentum reflection rates are asymmetric. In 
the mo — mjsi (m^ ~ mo) scattering, the negative (positive) mo- 
mentum reflection rates are obtained. The results explain the 
absence of thermal rectification in graphene: although the mo- 
mentum reflection rate is dependent upon the incident direc- 
tion, the energy reflection rate is directional independent. The 
same amount of energy would be reflected by the mass inter- 
face, thus it brings no reflection effect in the heat conduction 
process. 

Here we point out that the role of the KdV solitons sur- 
mised by Chang et al. plays in thermal rectification is still 
under heavy debate ifTl [T2m4l l37l . First, an extremely large 
excitation is needed to generate the supersonic KdV soli- 
tons. In carbon nanotubes, they are surmised to be generated 
by the collision of electrons, ultrashort laser light or strong 
compressions [38i .39] . So it is almost impossible that in the 
heat conduction process at room temperature, the KdV soli- 
tons could be generated. So far there is still no direct evi- 
dence of the supersonic KdV solitons in either carbon nan- 
otubes or graphene flakes. Second, Chang et al. suggest the 
prefeiTed direction of the heat flux is from the heavy to the 
light mass regions by considering KdV solitons lfTTI . How- 
ever in the MD simulations of carbon nanotubes, it shows that 
the preferred direction is from the light to the heavy mass 
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regions. It contradicts the existence of the KdV soHtons in 
carbon nanotube s lfT2l4T4ll . In our molecular dynamics simula- 
tions of graphene flakes, there is no thermal rectification effect 
in the heat conduction process. The result is also against the 
existence of the KdV solitons in graphene. Third, the square 
of the amplitudes is used to estimate the energy of a soliton 
by Chang et al. However, the dependence of energy upon the 
amplitude is far more complicated! 35, 36 1. Thus the width 
of a soliton is also necessary to be taken into consideration in 
order to measure the amount of energy. 

In summary, the inhomogeneous mass interface cannot 
bring thermal rectification effect in graphene. The absence 
of rectification effect is confirmed by studying different heat 
flux and mass difference ratio. The microscopic mechanism is 



explained by the interfacial scattering of solitons in graphene 
which reveals directional independent energy reflection rate. 
Our results imply that the mass interface or the mass gradient 
which is a combination of multiple mass interfaces cannot be 
applied to design graphene based thermal rectifiers. 
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